#setwd('/afs/inf.ed.ac.uk/user/s17/s1725186/Documents/PhD-Models/FirstPUModel/RMarkdowns')

library(tidyverse) ; library(reshape2) ; library(glue) ; library(plotly) ; library(plotlyutils)
library(RColorBrewer) ; library(viridis) ; require(gridExtra) ; library(GGally)
library(Rtsne)
library(ClusterR)
library(DESeq2)
library(expss)
library(knitr)

Load preprocessed dataset (preprocessing code in 19_10_14_data_preprocessing.Rmd)

# Gandal dataset
load('./../Data/preprocessed_data_imputed.RData')
datExpr = datExpr %>% data.frame
DE_info = DE_info %>% data.frame

# GO Neuronal annotations: regex 'neuron' in GO functional annotations and label the genes that make a match as neuronal
GO_annotations = read.csv('./../Data/genes_GO_annotations.csv')
GO_neuronal = GO_annotations %>% filter(grepl('neuron', go_term)) %>% 
              mutate('ID'=as.character(ensembl_gene_id)) %>% 
              dplyr::select(-ensembl_gene_id) %>% distinct(ID) %>%
              mutate('Neuronal'=1)

# SFARI Genes
SFARI_genes = read_csv('./../../../SFARI/Data/SFARI_genes_08-29-2019_with_ensembl_IDs.csv')
SFARI_genes = SFARI_genes[!duplicated(SFARI_genes$ID) & !is.na(SFARI_genes$ID),]

# Update DE_info with SFARI and Neuronal information
DE_info = DE_info %>% mutate('ID'=rownames(.)) %>% left_join(SFARI_genes, by='ID') %>% 
  mutate(`gene-score`=ifelse(is.na(`gene-score`), 'None', `gene-score`)) %>%
  distinct(ID, .keep_all = TRUE) %>% left_join(GO_neuronal, by='ID') %>%
  mutate(Neuronal=ifelse(is.na(Neuronal), 0, Neuronal)) %>%
  mutate(gene.score=ifelse(`gene-score`=='None' & Neuronal==1, 'Neuronal', `gene-score`), significant=padj<0.05 & !is.na(padj))


SFARI_colour_hue = function(r) {
  pal = c('#FF7631','#FFB100','#E8E328','#8CC83F','#62CCA6','#59B9C9','#b3b3b3','#808080','gray','#d9d9d9')[r]
}

SFARI Gene list

cat(paste0('There are ', length(unique(SFARI_genes$`gene-symbol`)), ' genes with a SFARI score'))
## There are 979 genes with a SFARI score

There are 979 genes with a SFARI score. but to map them to gene expression mapa we had to map the gene names to their corresponding ensembl IDs

Mapping SFARI Gene names to Ensembl IDs

cat(paste0('There are ', nrow(SFARI_genes), ' Ensembl IDs corresponding to the ',
             length(unique(SFARI_genes$`gene-symbol`)),' genes in the SFARI Gene dataset'))
## There are 1090 Ensembl IDs corresponding to the 979 genes in the SFARI Gene dataset



  • Since a gene can have more than one ensembl ID, there were some one-to-many mappings between a gene name and ensembl IDs, so that’s why we ended up with 1090 rows in the SFARI_genes dataset.

  • The details about how the genes were annotated with their Ensembl IDs can be found in 20_02_06_get_ensembl_ids.html

cat(paste0('There are ', sum(is.na(SFARI_genes$`gene-score`)) ,
             ' genes in the SFARI list without a score, of which ',
             sum(is.na(SFARI_genes$`gene-score`) & SFARI_genes$syndromic==0),
             ' don\'t have syndromic tag either (Why include them then???)'))
## There are 24 genes in the SFARI list without a score, of which 3 don't have syndromic tag either (Why include them then???)

Exploratory Analysis

cat(paste0('There are ', sum(SFARI_genes$ID %in% rownames(datExpr)), ' SFARI Genes in the expression dataset (~',
             round(100*mean(SFARI_genes$ID %in% rownames(datExpr))),'%)'))
## There are 924 SFARI Genes in the expression dataset (~85%)
cat(paste0('Of these, only ', sum(DE_info$`gene-score`!='None'), ' have an assigned score'))
## Of these, only 902 have an assigned score

From now on, we’re only going to focus on the 902 genes with a score

Gene count by SFARI score:

table_info = DE_info %>% apply_labels(`gene-score` = 'SFARI Gene Score', syndromic = 'Syndromic Tag',
                                      Neuronal = 'Neuronal Function', gene.score = 'Gene Score') %>%
             mutate(syndromic = as.logical(syndromic), Neuronal = as.logical(Neuronal))

cro(table_info$`gene-score`)
 #Total 
 SFARI Gene Score 
   1  25
   2  65
   3  191
   4  433
   5  165
   6  23
   None  15245
   #Total cases  16147

Gene count by Syndromic tag:

cro(table_info$syndromic)
 #Total 
 Syndromic Tag 
   FALSE  816
   TRUE  108
   #Total cases  924


GO Neuronal annotations:

cat(glue(sum(GO_neuronal$ID %in% rownames(datExpr)), ' genes have neuronal-related annotations'))
## 1094 genes have neuronal-related annotations
cat(glue(sum(SFARI_genes$ID %in% GO_neuronal$ID),' of these genes have a SFARI score'))
## 193 of these genes have a SFARI score
cro(table_info$gene.score[DE_info$`gene-score` %in% c('1','2','3','4','5','6')],
    list(table_info$Neuronal[DE_info$`gene-score` %in% c('1','2','3','4','5','6')], total()))
 Neuronal Function     #Total 
 FALSE   TRUE   
 Gene Score 
   1  16 9   25
   2  48 17   65
   3  149 42   191
   4  359 74   433
   5  129 36   165
   6  19 4   23
   #Total cases  720 182   902


Gene Expression


Normalised data

  • The higher the SFARI score, the higher the mean expression of the gene: This pattern is quite strong and it doesn’t have any biological interpretation, so it’s probably bias in the SFARI score assignment

  • The higher the SFARI score, the lower the standard deviation: This pattern is not as strong, but it is weird because the data was originally heteroscedastic with a positive relation between mean and variance, so the fact that the relation now seems to have reversed could mean that the vst normalisation ended up affecting the highly expressed genes more than it should have when trying to correct their higher variance

plot_data = data.frame('ID'=rownames(datExpr), 'MeanExpr'=rowMeans(datExpr), 'SDExpr'=apply(datExpr,1,sd)) %>% 
            left_join(DE_info, by='ID')

p1 = ggplotly(plot_data %>% ggplot(aes(gene.score, MeanExpr, fill=gene.score)) + geom_boxplot() + 
              scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + theme_minimal() +
              theme(legend.position='none'))

p2 = ggplotly(plot_data %>% ggplot(aes(gene.score, SDExpr, fill=gene.score)) + geom_boxplot() + 
              scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + theme_minimal() +
              ggtitle('Mean Expression (left) and SD (right) by SFARI score') + 
              theme(legend.position='none'))

subplot(p1, p2, nrows=1)
rm(plot_data, p1, p2)


Raw data

Just to corroborate that the relation between sd and SFARI score used to be in the opposite direction before the normalisation: The higher the SFARI score the higher the mean expression and the higher the standard deviation

*There are a lot of outliers, but the plot is interactive so you can zoom in

# Save preprocessed results
datExpr_prep = datExpr
datMeta_prep = datMeta
DE_info_prep = DE_info

load('./../Data/filtered_raw_data_imputed.RData')

plot_data = data.frame('ID'=rownames(datExpr), 'MeanExpr'=rowMeans(datExpr), 'SDExpr'=apply(datExpr,1,sd)) %>% 
            left_join(DE_info, by='ID')

p1 = ggplotly(plot_data %>% ggplot(aes(gene.score, MeanExpr, fill=gene.score)) + geom_boxplot() + 
              scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + theme_minimal() +
              theme(legend.position='none'))

p2 = ggplotly(plot_data %>% ggplot(aes(gene.score, SDExpr, fill=gene.score)) + geom_boxplot() + 
              scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + theme_minimal() +
              ggtitle('Mean Expression (left) and SD (right) by SFARI score') + 
              theme(legend.position='none'))

subplot(p1, p2, nrows=1)
rm(plot_data, p1, p2)

Return to normalised version of the data

# Save preprocessed results
datExpr = datExpr_prep
datMeta = datMeta_prep
DE_info = DE_info_prep

rm(datExpr_prep, datMeta_prep, DE_info_prep)


Log Fold Change

There seems to be a negative relation between SFARI score and log fold change when it would be expected to be either positively correlated or independent from each other (this last one because there are other factors that determine if a gene is releated to Autism apart from differences in gene expression)

Wikipedia mentions the likely explanation for this: “A disadvantage and serious risk of using fold change in this setting is that it is biased and may misclassify differentially expressed genes with large differences (B − A) but small ratios (B/A), leading to poor identification of changes at high expression levels”.

Based on this, since we saw there is a strong relation between SFARI score and mean expression, the bias in log fold change affects mainly genes with high SFARI scores, which would be the ones we are most interested in.

On top of this, I believe this effect is made more extreme by the pattern found in the previous plots, since the higher expressed genes were the most affected by the normalisation transformation, ending up with a smaller variance than the rest of the data, which is related to smaller ratios. (This is a constant problem independently of the normalisation function used).

ggplotly(DE_info %>% ggplot(aes(x=gene.score, y=abs(log2FoldChange), fill=gene.score)) + 
         geom_boxplot() + scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + 
         theme_minimal() + theme(legend.position='none'))


Effects of modifying filtering threshold by SFARI score

The higher the SFARI score, the higher the percentage of genes that get filtered by differential expression.

lfc_list = seq(1, 1.2, 0.01)

all_counts = data.frame('group'='All', 'n'=as.character(nrow(DE_info)))
Neuronal_counts = data.frame('group'='Neuronal', n=as.character(sum(DE_info$Neuronal)))
lfc_counts_all = DE_info %>% group_by(`gene-score`) %>% tally %>%
                 mutate('group'=as.factor(`gene-score`), 'n'=as.character(n)) %>%
                 dplyr::select(group, n) %>%
                 bind_rows(Neuronal_counts, all_counts) %>%
                 mutate('lfc'=-1) %>%  dplyr::select(lfc, group, n)

for(lfc in lfc_list){
  
  # Recalculate DE_info with the new threshold (p-values change)
  DE_genes = results(dds, lfcThreshold=log2(lfc), altHypothesis='greaterAbs') %>% data.frame
  
  DE_genes = DE_genes %>% mutate('ID'=rownames(.)) %>% left_join(SFARI_genes, by='ID') %>% 
             mutate(`gene-score`=ifelse(is.na(`gene-score`), 'None', `gene-score`)) %>%
             distinct(ID, .keep_all = TRUE) %>% left_join(GO_neuronal, by='ID') %>%
             mutate(Neuronal=ifelse(is.na(Neuronal), 0, Neuronal)) %>%
             mutate(gene.score=ifelse(`gene-score`=='None' & Neuronal==1, 'Neuronal', `gene-score`))
  
  DE_genes = DE_genes %>% filter(padj<0.05 & abs(log2FoldChange)>log2(lfc))

  
  # Calculate counts by groups
  all_counts = data.frame('group'='All', 'n'=as.character(nrow(DE_genes)))
  Neuronal_counts = data.frame('group'='Neuronal', n=as.character(sum(DE_genes$Neuronal)))
  lfc_counts = DE_genes %>% group_by(`gene-score`) %>% tally %>%
               mutate('group'=`gene-score`, 'n'=as.character(n)) %>%
               bind_rows(Neuronal_counts, all_counts) %>%
               mutate('lfc'=lfc) %>% dplyr::select(lfc, group, n)
  
  
  # Update lfc_counts_all
  lfc_counts_all = lfc_counts_all %>% bind_rows(lfc_counts)
}

# Add missing entries with 0s
lfc_counts_all = expand.grid('group'=unique(lfc_counts_all$group), 'lfc'=unique(lfc_counts_all$lfc)) %>% 
  left_join(lfc_counts_all, by=c('group','lfc')) %>% replace(is.na(.), 0)

# Calculate percentage of each group remaining
tot_counts = DE_info %>% group_by(`gene-score`) %>% tally() %>% filter(`gene-score`!='None') %>%
             mutate('group'=`gene-score`, 'tot'=n) %>% dplyr::select(group, tot) %>%
             bind_rows(data.frame('group'='Neuronal', 'tot'=sum(DE_info$Neuronal)),
                       data.frame('group'='All', 'tot'=nrow(DE_info)))

lfc_counts_all = lfc_counts_all %>% filter(lfc!=-1, group!='None') %>% 
                 left_join(tot_counts, by='group') %>% mutate('perc'=round(100*as.numeric(n)/tot,2))


# Plot change of number of genes
ggplotly(lfc_counts_all %>% ggplot(aes(lfc, perc, color=group)) + geom_point(aes(id=n)) + geom_line() + 
         scale_color_manual(values=SFARI_colour_hue(r=1:8)) + ylab('% of remaining genes') +  xlab('Fold Change') + 
         ggtitle('Effect of filtering thresholds by SFARI score') + theme_minimal())
rm(lfc_list, all_counts, Neuronal_counts, lfc_counts_all, lfc, lfc_counts, lfc_counts_all, tot_counts, lfc_counts_all)
cat(paste0('There are ', sum(DE_info$padj<0.05 & DE_info$`gene-score` != 'None' & !is.na(DE_info$padj)),
           ' SFARI genes that are differentially expressed'))
## There are 274 SFARI genes that are differentially expressed
kable(DE_info %>% filter(padj<0.05 & `gene-score` %in% c(1,2,3) & !is.na(padj)) %>% 
      dplyr::select(ID, `gene-symbol`, log2FoldChange, padj, `gene-score`, Neuronal) %>% arrange(`gene-score`,-abs(log2FoldChange)),
      caption = 'Top SFARI scores that are DE')
Top SFARI scores that are DE
ID gene-symbol log2FoldChange padj gene-score Neuronal
ENSG00000136535 TBR1 -0.2512818 0.0006964 1 1
ENSG00000197283 SYNGAP1 -0.1477068 0.0070142 1 1
ENSG00000136531 SCN2A -0.1351759 0.0052941 1 1
ENSG00000141431 ASXL3 0.1088269 0.0047548 1 0
ENSG00000110066 KMT5B 0.1038086 0.0001190 1 0
ENSG00000145362 ANK2 -0.0765338 0.0202311 1 1
ENSG00000174469 CNTNAP2 -0.2891387 0.0000087 2 1
ENSG00000169432 SCN9A 0.2814343 0.0102705 2 1
ENSG00000155974 GRIP1 -0.1976310 0.0017296 2 1
ENSG00000055609 KMT2C 0.1774044 0.0000296 2 0
ENSG00000114166 KAT2B 0.1423734 0.0016539 2 0
ENSG00000147050 KDM6A 0.1373929 0.0002086 2 0
ENSG00000061676 NCKAP1 -0.1333782 0.0000302 2 0
ENSG00000157103 SLC6A1 -0.1279723 0.0056354 2 1
ENSG00000108510 MED13 0.1247258 0.0000044 2 0
ENSG00000177030 DEAF1 -0.1236212 0.0018098 2 0
ENSG00000118482 PHF3 0.1205391 0.0012754 2 0
ENSG00000102786 INTS6 0.1189203 0.0064894 2 0
ENSG00000038382 TRIO 0.1070650 0.0004655 2 0
ENSG00000117139 KDM5B 0.1068038 0.0025796 2 0
ENSG00000169057 MECP2 -0.1031751 0.0173741 2 1
ENSG00000139613 SMARCC2 -0.0818166 0.0355780 2 0
ENSG00000187555 USP7 -0.0774175 0.0199852 2 0
ENSG00000141027 NCOR1 0.0755352 0.0162610 2 0
ENSG00000149930 TAOK2 -0.0654088 0.0292008 2 0
ENSG00000100354 TNRC6B 0.0559421 0.0457573 2 0
ENSG00000135046 ANXA1 0.6993722 0.0000503 3 0
ENSG00000171759 PAH -0.5479799 0.0001044 3 0
ENSG00000259207 ITGB3 0.4552313 0.0402947 3 0
ENSG00000144285 SCN1A -0.3866236 0.0000001 3 1
ENSG00000170745 KCNS3 -0.3198988 0.0007786 3 0
ENSG00000116117 PARD3B 0.2821225 0.0000003 3 0
ENSG00000157087 ATP2B2 -0.2766066 0.0000021 3 1
ENSG00000107099 DOCK8 0.2585677 0.0057421 3 0
ENSG00000196876 SCN8A -0.2583938 0.0000002 3 1
ENSG00000181722 ZBTB20 0.2571866 0.0000000 3 0
ENSG00000136854 STXBP1 -0.2511294 0.0000886 3 1
ENSG00000112655 PTK7 0.2439755 0.0068452 3 1
ENSG00000182621 PLCB1 -0.2392408 0.0000124 3 0
ENSG00000162946 DISC1 0.2359998 0.0002705 3 1
ENSG00000074590 NUAK1 -0.2342689 0.0000001 3 0
ENSG00000087085 ACHE -0.2260404 0.0080186 3 0
ENSG00000078328 RBFOX1 -0.2183744 0.0000034 3 0
ENSG00000166501 PRKCB -0.2165824 0.0005539 3 0
ENSG00000168769 TET2 0.2121650 0.0000004 3 0
ENSG00000124140 SLC12A5 -0.2004873 0.0006865 3 1
ENSG00000170579 DLGAP1 -0.1966584 0.0000340 3 0
ENSG00000197535 MYO5A -0.1924819 0.0000040 3 1
ENSG00000132294 EFR3A -0.1905926 0.0000006 3 0
ENSG00000128849 CGNL1 0.1821140 0.0156808 3 0
ENSG00000109911 ELP4 0.1811358 0.0006548 3 0
ENSG00000184156 KCNQ3 -0.1789998 0.0120325 3 0
ENSG00000183454 GRIN2A -0.1760026 0.0024942 3 1
ENSG00000181090 EHMT1 0.1723872 0.0000005 3 0
ENSG00000158321 AUTS2 0.1709908 0.0009373 3 0
ENSG00000205581 HMGN1 0.1709815 0.0008687 3 0
ENSG00000003147 ICA1 -0.1648788 0.0001663 3 0
ENSG00000176884 GRIN1 -0.1627296 0.0163288 3 1
ENSG00000175497 DPP10 -0.1626501 0.0002812 3 0
ENSG00000141646 SMAD4 0.1604987 0.0000931 3 1
ENSG00000152583 SPARCL1 -0.1572982 0.0199480 3 0
ENSG00000021645 NRXN3 -0.1551089 0.0010262 3 0
ENSG00000065526 SPEN 0.1533754 0.0099087 3 0
ENSG00000008083 JARID2 0.1384306 0.0056383 3 0
ENSG00000107077 KDM4C 0.1351744 0.0083689 3 0
ENSG00000149571 KIRREL3 0.1319958 0.0369379 3 1
ENSG00000166313 APBB1 -0.1292647 0.0002087 3 1
ENSG00000079482 OPHN1 0.1260252 0.0041773 3 0
ENSG00000005955 GGNBP2 -0.1248188 0.0000341 3 0
ENSG00000127914 AKAP9 0.1236977 0.0150098 3 0
ENSG00000152217 SETBP1 0.1185329 0.0322829 3 0
ENSG00000127124 HIVEP3 0.1161521 0.0444234 3 0
ENSG00000166833 NAV2 0.1158904 0.0331329 3 0
ENSG00000135439 AGAP2 -0.1149470 0.0380398 3 1
ENSG00000113328 CCNG1 0.1129688 0.0089413 3 1
ENSG00000146247 PHIP 0.1118427 0.0041948 3 0
ENSG00000106290 TAF6 -0.1111955 0.0304439 3 0
ENSG00000197102 DYNC1H1 -0.1110523 0.0022439 3 0
ENSG00000104388 RAB2A -0.1076087 0.0003218 3 0
ENSG00000185345 PRKN -0.1074080 0.0159131 3 1
ENSG00000132604 TERF2 -0.1057978 0.0033537 3 0
ENSG00000142230 SAE1 -0.1007365 0.0026293 3 0
ENSG00000083168 KAT6A 0.1002845 0.0005914 3 0
ENSG00000092964 DPYSL2 -0.0991571 0.0066131 3 1
ENSG00000117362 APH1A 0.0990303 0.0475002 3 0
ENSG00000139726 DENR -0.0977187 0.0014690 3 0
ENSG00000164506 STXBP5 -0.0975439 0.0191122 3 0
ENSG00000095564 BTAF1 0.0936371 0.0166415 3 0
ENSG00000103197 TSC2 -0.0830945 0.0027326 3 1
ENSG00000146830 GIGYF1 -0.0821281 0.0405157 3 0
ENSG00000050030 NEXMIF -0.0791094 0.0468104 3 0
ENSG00000151150 ANK3 -0.0782771 0.0410641 3 1
ENSG00000182771 GRID1 -0.0715553 0.0464685 3 0

Session info

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] knitr_1.28                  expss_0.10.2               
##  [3] DESeq2_1.24.0               SummarizedExperiment_1.14.1
##  [5] DelayedArray_0.10.0         BiocParallel_1.18.1        
##  [7] matrixStats_0.56.0          Biobase_2.44.0             
##  [9] GenomicRanges_1.36.1        GenomeInfoDb_1.20.0        
## [11] IRanges_2.18.3              S4Vectors_0.22.1           
## [13] BiocGenerics_0.30.0         ClusterR_1.2.1             
## [15] gtools_3.8.2                Rtsne_0.15                 
## [17] GGally_1.5.0                gridExtra_2.3              
## [19] viridis_0.5.1               viridisLite_0.3.0          
## [21] RColorBrewer_1.1-2          plotlyutils_0.0.0.9000     
## [23] plotly_4.9.2                glue_1.3.2                 
## [25] reshape2_1.4.3              forcats_0.5.0              
## [27] stringr_1.4.0               dplyr_0.8.5                
## [29] purrr_0.3.3                 readr_1.3.1                
## [31] tidyr_1.0.2                 tibble_3.0.0               
## [33] ggplot2_3.3.0               tidyverse_1.3.0            
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.4-1       ellipsis_0.3.0         htmlTable_1.13.3      
##  [4] XVector_0.24.0         base64enc_0.1-3        fs_1.4.0              
##  [7] rstudioapi_0.11        bit64_0.9-7            AnnotationDbi_1.46.1  
## [10] fansi_0.4.1            lubridate_1.7.4        xml2_1.2.5            
## [13] splines_3.6.3          geneplotter_1.62.0     Formula_1.2-3         
## [16] jsonlite_1.6.1         annotate_1.62.0        broom_0.5.5           
## [19] cluster_2.1.0          dbplyr_1.4.2           png_0.1-7             
## [22] compiler_3.6.3         httr_1.4.1             backports_1.1.5       
## [25] assertthat_0.2.1       Matrix_1.2-18          lazyeval_0.2.2        
## [28] cli_2.0.2              acepack_1.4.1          htmltools_0.4.0       
## [31] tools_3.6.3            gmp_0.5-13.6           gtable_0.3.0          
## [34] GenomeInfoDbData_1.2.1 Rcpp_1.0.4             cellranger_1.1.0      
## [37] vctrs_0.2.4            nlme_3.1-144           crosstalk_1.1.0.1     
## [40] xfun_0.12              rvest_0.3.5            lifecycle_0.2.0       
## [43] XML_3.99-0.3           zlibbioc_1.30.0        scales_1.1.0          
## [46] hms_0.5.3              yaml_2.2.1             memoise_1.1.0         
## [49] rpart_4.1-15           RSQLite_2.2.0          reshape_0.8.8         
## [52] latticeExtra_0.6-29    stringi_1.4.6          highr_0.8             
## [55] genefilter_1.66.0      checkmate_2.0.0        rlang_0.4.5           
## [58] pkgconfig_2.0.3        bitops_1.0-6           evaluate_0.14         
## [61] lattice_0.20-40        labeling_0.3           htmlwidgets_1.5.1     
## [64] bit_1.1-15.2           tidyselect_1.0.0       plyr_1.8.6            
## [67] magrittr_1.5           R6_2.4.1               generics_0.0.2        
## [70] Hmisc_4.4-0            DBI_1.1.0              pillar_1.4.3          
## [73] haven_2.2.0            foreign_0.8-75         withr_2.1.2           
## [76] survival_3.1-11        RCurl_1.98-1.1         nnet_7.3-13           
## [79] modelr_0.1.6           crayon_1.3.4           rmarkdown_2.1         
## [82] jpeg_0.1-8.1           locfit_1.5-9.4         grid_3.6.3            
## [85] readxl_1.3.1           data.table_1.12.8      blob_1.2.1            
## [88] reprex_0.3.0           digest_0.6.25          xtable_1.8-4          
## [91] munsell_0.5.0